Data Preparation

packages = c('tmap', 'SpatialAcc', 'sf', 'ggstatsplot', 'tidyverse')

for(p in packages){
  if(!require(p, character.only = T)){
    install.packages(p)
  }
  library(p, character.only = T)
}

1 Importing Geospatial Data

na_validity_check <-function(input){
  print("Check for CRS...")
  print(st_crs(input))
  print("Check for NA...")
  print(input[rowSums(is.na(input))!=0,])
  test <- st_is_valid(input,reason=TRUE)
  print("Check for invalid geometry...")
  length(which(test!= "Valid Geometry"))
}

subzone = st_read(dsn = "data/geospatial/subzone.kml")

## Reading layer `URA_MP19_SUBZONE_NO_SEA_PL' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\subzone.kml' using driver `KML'
## Simple feature collection with 332 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XYZ
## bbox:           xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

cib = st_read(dsn = "data/geospatial", layer = "cib_gardens")

## Reading layer `cib_gardens' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 1024 features and 14 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 5781.933 ymin: 27790.41 xmax: 45160.86 ymax: 48982.52
## projected CRS:  SVY21

comm_club = st_read(dsn = "data/geospatial/community_clubs.kml")

## Reading layer `COMMUNITYCLUBS' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\community_clubs.kml' using driver `KML'
## Simple feature collection with 119 features and 2 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: 103.6923 ymin: 1.274863 xmax: 103.9592 ymax: 1.451682
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

comm_use_site = st_read(dsn = "data/geospatial/community_use_sites.kml")

## Reading layer `COMM_USE_SITES' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\community_use_sites.kml' using driver `KML'
## Simple feature collection with 250 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XYZ
## bbox:           xmin: 103.6228 ymin: 1.269832 xmax: 103.9653 ymax: 1.461749
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

dus_sports_facilities = st_read(dsn = "data/geospatial", layer = "dus_school_sports_facilities")

## Reading layer `dus_school_sports_facilities' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 174 features and 15 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 11739.27 ymin: 28604.3 xmax: 42410.51 ymax: 48696.35
## projected CRS:  SVY21

nature_areas = st_read(dsn = "data/geospatial/nature_areas.kml")

## Reading layer `UP_G_MP19_PKWB_NATURE_PL' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\nature_areas.kml' using driver `KML'
## Simple feature collection with 62 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XYZ
## bbox:           xmin: 103.644 ymin: 1.247164 xmax: 104.0767 ymax: 1.451329
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

np_activity_area = st_read(dsn = "data/geospatial/nparks_activity_area.kml")

## Reading layer `NPARKS_ACTIVITY_AREA' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\nparks_activity_area.kml' using driver `KML'
## Simple feature collection with 599 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XYZ
## bbox:           xmin: 103.693 ymin: 1.262527 xmax: 103.9936 ymax: 1.462811
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

np_parks = st_read(dsn = "data/geospatial/nparks_parks.kml")

## Reading layer `NPARKS_PARKS' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\nparks_parks.kml' using driver `KML'
## Simple feature collection with 390 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XYZ
## bbox:           xmin: 103.6925 ymin: 1.265889 xmax: 104.008 ymax: 1.46419
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

np_play_fitness = st_read(dsn = "data/geospatial/nparks_play_fitness_equipment.kml")

## Reading layer `NPARKS_PLAY_FITNESS_EQ' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\nparks_play_fitness_equipment.kml' using driver `KML'
## Simple feature collection with 3108 features and 2 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: 103.6929 ymin: 1.262622 xmax: 103.9935 ymax: 1.457533
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

park_loop = st_read(dsn = "data/geospatial/park_connector_loop.kml")

## Reading layer `PARK_CONNECTOR_LOOP' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\park_connector_loop.kml' using driver `KML'
## Simple feature collection with 233 features and 2 fields
## geometry type:  LINESTRING
## dimension:      XYZ
## bbox:           xmin: 103.6933 ymin: 1.272949 xmax: 104.0316 ymax: 1.462852
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

preschools = st_read(dsn = "data/geospatial/preschools.kml")

## Reading layer `PRESCHOOLS_LOCATION' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\preschools.kml' using driver `KML'
## Simple feature collection with 1925 features and 2 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: 103.6824 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

schools = st_read(dsn = "data/geospatial", layer = "schools")

## Reading layer `schools' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 327 features and 43 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 11799.7 ymin: 28571.39 xmax: 42434.9 ymax: 48727.31
## projected CRS:  SVY21 / Singapore TM

student_care = st_read(dsn = "data/geospatial/student_care.kml")

## Reading layer `STUDENTCARE' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\student_care.kml' using driver `KML'
## Simple feature collection with 425 features and 2 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: 103.6877 ymin: 1.27474 xmax: 103.963 ymax: 1.456667
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

sportsg_facilities= st_read(dsn = "data/geospatial/sportsg_sports_facilities.kml")

## Reading layer `SPORTSG_SPORT_FACILITIES' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\sportsg_sports_facilities.kml' using driver `KML'
## Simple feature collection with 35 features and 2 fields
## geometry type:  POLYGON
## dimension:      XYZ
## bbox:           xmin: 103.6932 ymin: 1.285977 xmax: 103.9526 ymax: 1.435855
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

watersports_facilities = st_read(dsn = "data/geospatial", layer = "water_sport_facilities")

## Reading layer `water_sport_facilities' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 17 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 12469.6 ymin: 30191.07 xmax: 45423.18 ymax: 48969.16
## projected CRS:  SVY21

1.1 Data Cleaning and Transforming Subzone Layer

na_validity_check(subzone)

## [1] "Check for CRS..."
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
## [1] "Check for NA..."
## Simple feature collection with 0 features and 2 fields
## bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
## geographic CRS: WGS 84
## [1] Name        Description geometry   
## <0 rows> (or 0-length row.names)
## [1] "Check for invalid geometry..."

## [1] 6

subzone <- st_set_crs(subzone, 3414)
subzone<- st_make_valid(subzone)

getName <- function(input){
  noHtml <- str_replace_all(input,"<.*?>","")
  matched_str <- str_match(noHtml, "\\sSUBZONE_N (.*?) SUBZONE_C")[,2]
  return(str_trim(matched_str))
}

subzone <- subzone %>%
  mutate(`SUBZONE_N` = getName(subzone$Description)) 
head(subzone, 1)

1.2 Data transformation for required layers

nature_areas <- st_transform(nature_areas, 3414)
nature_areas<- st_make_valid(nature_areas)
np_activity_area <- st_transform(np_activity_area, 3414)
np_parks <- st_transform(np_parks, 3414)
np_parks<- st_make_valid(np_parks)
np_play_fitness <- st_transform(np_play_fitness, 3414)
park_loop <- st_transform(park_loop, 3414)
preschools <- st_transform(preschools, 3414)
sportsg_facilities <- st_transform(sportsg_facilities, 3414)
student_care <- st_transform(student_care, 3414)
watersports_facilities <- st_set_crs(watersports_facilities, 3414)
cib <- st_set_crs(cib, 3414)
comm_club <- st_transform(comm_club, 3414)
comm_use_site <- st_transform(comm_use_site, 3414)
dus_sports_facilities <- st_set_crs(dus_sports_facilities, 3414)

2 Import Aspatial Data

HDB

hdb <- read_csv('data/aspatial/hdb_osm.csv')

## Parsed with column specification:
## cols(
##   .default = col_logical(),
##   X = col_double(),
##   Y = col_double(),
##   full_id = col_character(),
##   osm_id = col_double(),
##   osm_type = col_character(),
##   addr_city = col_character(),
##   addr_house = col_character(),
##   addr_postc = col_character(),
##   addr_stree = col_character(),
##   building = col_character(),
##   building_l = col_character(),
##   name = col_character(),
##   type = col_character(),
##   addr_count = col_character(),
##   residentia = col_character(),
##   year = col_double(),
##   addr_hou_1 = col_character(),
##   building_c = col_character(),
##   building_u = col_character(),
##   building_m = col_character()
##   # ... with 41 more columns
## )

## See spec(...) for full column specifications.

## Warning: 819 parsing failures.
##  row        col           expected       actual                        file
## 1029 amenity    1/0/T/F/TRUE/FALSE parking      'data/aspatial/hdb_osm.csv'
## 1029 parking    1/0/T/F/TRUE/FALSE multi-storey 'data/aspatial/hdb_osm.csv'
## 1193 operator_2 1/0/T/F/TRUE/FALSE HDB          'data/aspatial/hdb_osm.csv'
## 1194 operator_2 1/0/T/F/TRUE/FALSE HDB          'data/aspatial/hdb_osm.csv'
## 1195 operator_2 1/0/T/F/TRUE/FALSE HDB          'data/aspatial/hdb_osm.csv'
## .... .......... .................. ............ ...........................
## See problems(...) for more details.

hdb <- st_as_sf(hdb, 
                coords = c('X', 'Y'),
                crs = 4326) %>%
  st_transform(3414)

population <- read_csv('data/aspatial/population.csv')

## Parsed with column specification:
## cols(
##   PA = col_character(),
##   SZ = col_character(),
##   AG = col_character(),
##   Sex = col_character(),
##   TOD = col_character(),
##   Pop = col_double(),
##   Time = col_double()
## )

3 Data preparation

• Extract population data for 2020 only
• Calculate children population living in HDB: children aged 0 to 7
  • Number of children aged 5 to 7 is estimated as 3/5 of the total number of children aged 5 to 9

\[Number\ of\ children\ aged\ 0\ to\ 7\ in\ subzone\ =\ Number\ of\ children\ aged\ 0\ to\ 4 + (\frac{3}{5}*Number\ of\ children\ aged\ 5\ to\ 9)\]

children_population <- population %>%
  filter(Time == 2020) %>%
  group_by(PA, SZ, AG, TOD) %>%
  summarise(`POP` = sum(`Pop`)) %>%
  ungroup() %>%
  pivot_wider(names_from = AG, values_from = POP) %>%
  mutate(children_population = `0_to_4` + (3/5 * `5_to_9`)) %>%
  select(PA, SZ, TOD, children_population) %>%
  pivot_wider(names_from = TOD, values_from = children_population) %>%
  mutate_at(.vars = vars(PA, SZ), .funs = ~ toupper(.)) %>%
  mutate(children_hdb = rowSums(.[4:7])) %>%
  select(SZ, children_hdb)

## `summarise()` regrouping output by 'PA', 'SZ', 'AG' (override with `.groups` argument)

children_population

Join population data with subzone data

subzone_children <- subzone %>%
  left_join(children_population, by = c('SUBZONE_N' = 'SZ')) %>%
  select(-c(Description, Name))

glimpse(subzone_children)

## Rows: 332
## Columns: 3
## $ SUBZONE_N    <chr> "MARINA EAST", "INSTITUTION HILL", "ROBERTSON QUAY", "...
## $ children_hdb <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 492, 0, 358,...
## $ geometry     <GEOMETRY [m]> MULTIPOLYGON (((103.8802 1...., MULTIPOLYGON ...

Calculate number of children aged 0 to 7 in each HDB block

• Children aged 0 to 7 are the key focus of this project, as these are the formative years of children development.
• The number of children residing in each HDB block of a subzone will be estimated by dividing the total number of children living in the subzone, by the total number of HDB blocks in that subzone.

\[Number\ of\ children\ in\ each\ HDB\ =\ \frac{Number\ of\ children\ in\ subzone}{Number\ of\ HDB\ blocks\ in\ subzone}\]

Count number of HDB blocks in each subzone

subzone_children$num_hdb <- lengths(st_intersects(subzone_children, hdb))

summary(subzone_children$num_hdb)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0

qtm(subzone_children, 'num_hdb')

qtm(subzone_children, 'children_hdb')

Calculate number of children in each HDB: demand

subzone_children <- subzone_children %>%
  # Calculate demand. If there are no HDB
  mutate(children_each_hdb = ifelse(num_hdb != 0, children_hdb / num_hdb, 0))

summary(subzone_children$children_each_hdb)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0

tm_shape(subzone_children) +
  tm_polygons('children_each_hdb',
              style = 'jenks')

Add information on number of children for each HDB to HDB data

hdb <- hdb %>%
  select(`ID`, `addr_house`, `addr_postc`, `addr_stree`, `building_l`) %>%
  rename(house_number = addr_house,
         postal_code = addr_postc,
         street = addr_stree,
         num_levels = building_l) %>%
  st_join(subzone_children, join = st_intersects)

glimpse(hdb)

## Rows: 9,806
## Columns: 10
## $ ID                <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15...
## $ house_number      <chr> "21", "21", "94B", "94B", "94C", "94C", "25B", "2...
## $ postal_code       <chr> "270021", "270021", "461094", "461094", "462094",...
## $ street            <chr> "Ghim Moh Road", "Ghim Moh Road", "Bedok North Av...
## $ num_levels        <chr> "10", "10", "13", "13", "13", "13", "27", "27", N...
## $ SUBZONE_N         <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N...
## $ children_hdb      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N...
## $ num_hdb           <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N...
## $ children_each_hdb <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N...
## $ geometry          <POINT [m]> POINT (23015.11 32522.95), POINT (23015.11 ...

4 Calculate accessibility

4.1 Student care

student_care <- st_read(dsn = 'data/geospatial/distance_matrix_', layer = 'student_care') %>%
  select(ID, Name, ADDRESSSTR, ADDRESSPOS) %>%
  rename(name = Name,
         address = ADDRESSSTR,
         postal_code = ADDRESSPOS)

## Reading layer `student_care' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\distance_matrix_' using driver `ESRI Shapefile'
## Simple feature collection with 425 features and 28 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 11796.89 ymin: 28579.85 xmax: 42433.39 ymax: 48696.35
## projected CRS:  SVY21 / Singapore TM

glimpse(student_care)

## Rows: 425
## Columns: 5
## $ ID          <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ...
## $ name        <chr> "Ananias Centre (Choa Chu Kang)", "Ananias Centre (Clem...
## $ address     <chr> "Blk 289 Choa Chu Kang Avenue 3 #01-262 SINGAPORE 68028...
## $ postal_code <chr> "680289", "121440", "760160", "760417", "530540", "6801...
## $ geometry    <POINT [m]> POINT (17694.75 40144.54), POINT (20305.83 33142....

dm_student_care <- read_csv('data/aspatial/distance matrix/hdb_studentcare.csv')

## Parsed with column specification:
## cols(
##   origin_id = col_double(),
##   destination_id = col_double(),
##   entry_cost = col_double(),
##   network_cost = col_double(),
##   exit_cost = col_double(),
##   total_cost = col_double()
## )

glimpse(dm_student_care)

## Rows: 4,167,550
## Columns: 6
## $ origin_id      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ destination_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ entry_cost     <dbl> 89.29121, 89.29121, 89.29121, 89.29121, 89.29121, 89...
## $ network_cost   <dbl> 13126.843, 13126.843, 19604.224, 19672.647, 17050.73...
## $ exit_cost      <dbl> 35.10268, 35.10268, 36.57691, 98.50361, 146.68192, 3...
## $ total_cost     <dbl> 13251.237, 13251.237, 19730.092, 19860.441, 17286.70...

Tidy matrix – make wider

dm_student_care <- dm_student_care %>%
  # Obtain final matrix
  select(origin_id, destination_id, total_cost) %>%
  pivot_wider(names_from = destination_id, values_from = total_cost) %>%
  select(-origin_id)
  
# Convert to km
dm_student_care <- as.matrix(dm_student_care/1000)
# dm_student_care

Calculate average student care capacity for each student care centre

• Total capacity of all student care centres in Singapore is 42,907 in 2019.
• Total number of student care centres in Singapore is 425.

\[Average\ student\ care\ capacity\ for\ each\ centre\ = \frac{Total\ capacity\ of\ student\ care\ in\ Singapore}{Number\ of\ student\ care\ centres\ in\ Singapore}\]

student_care_capacity <- round(42907 / nrow(student_care), 0)
student_care_capacity

## [1] 101

• Average student care capacity for each centre is estimated to be 101.

Calculate Hansen accessibility

student_care_capacity_list <- rep(student_care_capacity, nrow(student_care))

acc_Hansen <- data.frame(ac(hdb$children_each_hdb,
                  student_care_capacity_list,
                  dm_student_care, 
                  d0 = 50,
                  power = 1.5, 
                  family = "Hansen"))

colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hansen_student_care <- bind_cols(hdb, acc_Hansen)

hansen_student_care

4.2 School

Import schools

pri_schools <- st_read(dsn = 'data/geospatial/distance_matrix_', layer = 'schools_primary') %>%
  select(ID, school_nam, address, postal_cod) %>%
  rename(name = school_nam,
         postal_code = postal_cod)

## Reading layer `schools_primary' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\distance_matrix_' using driver `ESRI Shapefile'
## Simple feature collection with 170 features and 44 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 11799.7 ymin: 28571.39 xmax: 42434.9 ymax: 48727.31
## projected CRS:  SVY21 / Singapore TM

glimpse(pri_schools)

## Rows: 170
## Columns: 5
## $ ID          <dbl> 1, 3, 5, 6, 7, 10, 16, 17, 18, 22, 24, 28, 30, 31, 37, ...
## $ name        <chr> "ADMIRALTY PRIMARY SCHOOL", "AHMAD IBRAHIM PRIMARY SCHO...
## $ address     <chr> "11   WOODLANDS CIRCLE", "10   YISHUN STREET 11", "100 ...
## $ postal_code <chr> "738907", "768643", "579646", "159016", "544969", "5699...
## $ geometry    <POINT [m]> POINT (24330.56 47178.67), POINT (27932.25 46173....

Import distance matrix

dm_pri_schools <- read_csv('data/aspatial/distance matrix/hdb_school.csv')

## Parsed with column specification:
## cols(
##   origin_id = col_double(),
##   destination_id = col_double(),
##   entry_cost = col_double(),
##   network_cost = col_double(),
##   exit_cost = col_double(),
##   total_cost = col_double()
## )

glimpse(dm_pri_schools)

## Rows: 1,667,020
## Columns: 6
## $ origin_id      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ destination_id <dbl> 1, 3, 5, 6, 7, 10, 16, 17, 18, 22, 24, 28, 30, 31, 3...
## $ entry_cost     <dbl> 89.29121, 89.29121, 89.29121, 89.29121, 89.29121, 89...
## $ network_cost   <dbl> 19430.115, 19880.538, 11378.661, 6226.051, 18680.708...
## $ exit_cost      <dbl> 88.249817, 56.515963, 80.705633, 7.244451, 93.627581...
## $ total_cost     <dbl> 19607.656, 20026.345, 11548.657, 6322.587, 18863.627...

Tidy matrix – make wider

dm_pri_schools <- dm_pri_schools %>%
  # Obtain final matrix
  select(origin_id, destination_id, total_cost) %>%
  pivot_wider(names_from = destination_id, values_from = total_cost) %>%
  select(-origin_id)
  
# Convert to km
dm_pri_schools <- as.matrix(dm_pri_schools/1000)

Calculate average primary school capacity for each school

• Enrolment for primary 1 will be utilised to calculate primary school capacity.
• Total number of primary schools in Singapore is 170.
• Primary 1 enrolment for 2019 is 37,671

school_capacity <- round(37671 / nrow(pri_schools))
school_capacity

## [1] 222

Calculate Hansen accessibility

school_capacity_list <- rep(school_capacity, nrow(pri_schools))

acc_Hansen <- data.frame(ac(hdb$children_each_hdb,
                  school_capacity_list,
                  dm_pri_schools, 
                  d0 = 50,
                  power = 2, 
                  family = "Hansen"))

colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hansen_pri_schools <- bind_cols(hdb, acc_Hansen)

hansen_pri_schools

4.3 Water sports facilities

Import water sports facilities

watersports_facilities <- st_read(dsn = 'data/geospatial/distance_matrix_', layer = 'water_sport_facilities') %>%
  select(ID, NAME, ADDRESSSTR, ADDRESSPOS) %>%
  rename(name = NAME,
         address = ADDRESSSTR,
         postal_code = ADDRESSPOS)

## Reading layer `water_sport_facilities' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\distance_matrix_' using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 18 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 12469.6 ymin: 30191.07 xmax: 45423.18 ymax: 48969.16
## projected CRS:  SVY21 / Singapore TM

glimpse(watersports_facilities)

## Rows: 32
## Columns: 5
## $ ID          <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ...
## $ name        <chr> "Toa Payoh Sports and Recreation Centre", "Water-Ventur...
## $ address     <chr> "Toa Payoh Lor. 6", "East Coast Parkway", "Bishan Stree...
## $ postal_code <int> 319392, 468961, 579778, 648965, 339627, 436752, 148948,...
## $ geometry    <POINT [m]> POINT (29869.2 34746.17), POINT (41233.04 32670.0...

Import distance matrix

dm_watersports_facilities <- read_csv('data/aspatial/distance matrix/hdb_watersports.csv')

## Parsed with column specification:
## cols(
##   origin_id = col_double(),
##   destination_id = col_double(),
##   entry_cost = col_double(),
##   network_cost = col_double(),
##   exit_cost = col_double(),
##   total_cost = col_double()
## )

glimpse(dm_watersports_facilities)

## Rows: 313,792
## Columns: 6
## $ origin_id      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ destination_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ entry_cost     <dbl> 89.29121, 89.29121, 89.29121, 89.29121, 89.29121, 89...
## $ network_cost   <dbl> 9060.498, 20286.510, 10517.724, 10517.724, 11776.135...
## $ exit_cost      <dbl> 57.776199, 38.678467, 86.974437, 86.974437, 31.99397...
## $ total_cost     <dbl> 9207.565, 20414.480, 10693.990, 10693.990, 11897.420...

Tidy matrix – make wider

dm_watersports_facilities <- dm_watersports_facilities %>%
  # Obtain final matrix
  select(origin_id, destination_id, total_cost) %>%
  pivot_wider(names_from = destination_id, values_from = total_cost) %>%
  select(-origin_id)
  
# Convert to km
dm_watersports_facilities <- as.matrix(dm_watersports_facilities/1000)

Calculate Hansen accessibility

watersports_capacity_list <- rep(1, nrow(watersports_facilities))

acc_Hansen <- data.frame(ac(hdb$children_each_hdb,
                  watersports_capacity_list,
                  dm_watersports_facilities, 
                  d0 = 50,
                  power = 1.5, 
                  family = "Hansen"))

colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hansen_watersports_facilities <- bind_cols(hdb, acc_Hansen)

hansen_watersports_facilities

5 Dual-use scheme school sports facilities

Import DUS sports facilities

dus_sports_facilities <- st_read(dsn = 'data/geospatial/distance_matrix_', layer = 'dus_school_sports_facilities') %>%
  select(ID, SCHOOL_NAM, ADDRESS, POSTAL_COD, FACILITIES) %>%
  rename(name = SCHOOL_NAM,
         address = ADDRESS,
         postal_code = POSTAL_COD,
         facilities = FACILITIES)

## Reading layer `dus_school_sports_facilities' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\distance_matrix_' using driver `ESRI Shapefile'
## Simple feature collection with 174 features and 16 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 11739.27 ymin: 28604.3 xmax: 42410.51 ymax: 48696.35
## projected CRS:  SVY21 / Singapore TM

glimpse(dus_sports_facilities)

## Rows: 174
## Columns: 6
## $ ID          <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ...
## $ name        <chr> "Changkat Changi Secondary School", "Chestnut Drive Sec...
## $ address     <chr> "23 Simei Street 3", "58 Chestnut Drive", "20 Choa Chu ...
## $ postal_code <chr> "529894", "679301", "689905", "688845", "129903", "1299...
## $ facilities  <chr> "Chargeable Field", "Chargeable Field, Basketball Court...
## $ geometry    <POINT [m]> POINT (41288.67 35836.16), POINT (21042.42 38992....

Import distance matrix

dm_dus_school_sports_facilities <- read_csv('data/aspatial/distance matrix/hdb_dus.csv')

## Parsed with column specification:
## cols(
##   origin_id = col_double(),
##   destination_id = col_double(),
##   entry_cost = col_double(),
##   network_cost = col_double(),
##   exit_cost = col_double(),
##   total_cost = col_double()
## )

glimpse(dm_dus_school_sports_facilities)

## Rows: 1,706,244
## Columns: 6
## $ origin_id      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ destination_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ entry_cost     <dbl> 89.29121, 89.29121, 89.29121, 89.29121, 89.29121, 89...
## $ network_cost   <dbl> 22496.583, 9295.449, 13681.212, 11848.690, 11848.690...
## $ exit_cost      <dbl> 100.22417, 45.43876, 80.20923, 39.46815, 39.46815, 1...
## $ total_cost     <dbl> 22686.099, 9430.179, 13850.712, 11977.449, 11977.449...

Tidy matrix – make wider

dm_dus_school_sports_facilities <- dm_dus_school_sports_facilities %>%
  # Obtain final matrix
  select(origin_id, destination_id, total_cost) %>%
  pivot_wider(names_from = destination_id, values_from = total_cost) %>%
  select(-origin_id)
  
# Convert to km
dm_dus_school_sports_facilities <- as.matrix(dm_dus_school_sports_facilities/1000)

Calculate Hansen accessibility

dus_capacity_list <- rep(1, nrow(dus_sports_facilities))

acc_Hansen <- data.frame(ac(hdb$children_each_hdb,
                  dus_capacity_list,
                  dm_dus_school_sports_facilities, 
                  d0 = 50,
                  power = 2, 
                  family = "Hansen"))

colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hansen_dus_school_sports_facilities <- bind_cols(hdb, acc_Hansen)

hansen_dus_school_sports_facilities

6 Community in bloom gardens

Import community in bloom gardens

cib_gardens <- st_read(dsn = 'data/geospatial/distance_matrix_', layer = 'cib_gardens') %>%
  select(ID, ADDRESS, DIVISION, CONSTITUEN, CATEGORY) %>%
  rename(address = ADDRESS,
         division = DIVISION, 
         constituency = CONSTITUEN,
         category = CATEGORY)

## Reading layer `cib_gardens' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\distance_matrix_' using driver `ESRI Shapefile'
## Simple feature collection with 1024 features and 15 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 5781.933 ymin: 27790.41 xmax: 45160.86 ymax: 48982.52
## projected CRS:  SVY21 / Singapore TM

glimpse(cib_gardens)

## Rows: 1,024
## Columns: 6
## $ ID           <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,...
## $ address      <chr> "Blk 730 Tampines St 71 #01-61", "TAI HWAN TERRACE, BL...
## $ division     <chr> "TAMPINES CENTRAL", "SERANGOON", "TAMPINES CENTRAL", "...
## $ constituency <chr> "TAMPINES GRC", "ALJUNIED GRC", "TAMPINES GRC", "TAMPI...
## $ category     <chr> "Public Estate", "Private Estate", "Public Estate", "P...
## $ geometry     <POINT [m]> POINT (39151.3 37639.3), POINT (31058.07 37624.8...

Import distance matrix

dm_cib_gardens <- read_csv('data/aspatial/distance matrix/hdb_cib.csv')

## Parsed with column specification:
## cols(
##   origin_id = col_double(),
##   destination_id = col_double(),
##   entry_cost = col_double(),
##   network_cost = col_double(),
##   exit_cost = col_double(),
##   total_cost = col_double()
## )

glimpse(dm_cib_gardens)

## Rows: 10,041,344
## Columns: 6
## $ origin_id      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ destination_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ entry_cost     <dbl> 89.29121, 89.29121, 89.29121, 89.29121, 89.29121, 89...
## $ network_cost   <dbl> 19708.18, 13155.52, 19763.04, 21985.90, 21768.26, 22...
## $ exit_cost      <dbl> 44.587488, 4.749120, 109.933316, 16.257711, 67.14973...
## $ total_cost     <dbl> 19842.06, 13249.56, 19962.26, 22091.45, 21924.70, 22...

Tidy matrix – make wider

dm_cib_gardens <- dm_cib_gardens %>%
  # Obtain final matrix
  select(origin_id, destination_id, total_cost) %>%
  pivot_wider(names_from = destination_id, values_from = total_cost) %>%
  select(-origin_id)
  
# Convert to km
dm_cib_gardens <- as.matrix(dm_cib_gardens/1000)

Calculate Hansen accessibility

cib_capacity_list <- rep(1, nrow(cib_gardens))

acc_Hansen <- data.frame(ac(hdb$children_each_hdb,
                  cib_capacity_list,
                  dm_cib_gardens, 
                  d0 = 50,
                  power = 2.5, 
                  family = "Hansen"))

colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hansen_cib_gardens <- bind_cols(hdb, acc_Hansen)

hansen_cib_gardens

7 Pre-schools

Import pre-schools

preschools <- st_read(dsn = 'data/geospatial/distance_matrix_', layer = 'preschools') %>%
  select(ID, CENTRE_NAM, ADDRESS, POSTAL_COD) %>%
  rename(name = CENTRE_NAM,
         address = ADDRESS,
         postal_code = POSTAL_COD)

## Reading layer `preschools' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\distance_matrix_' using driver `ESRI Shapefile'
## Simple feature collection with 1925 features and 19 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 11203.01 ymin: 25596.33 xmax: 45404.24 ymax: 49300.88
## projected CRS:  SVY21 / Singapore TM

glimpse(preschools)

## Rows: 1,925
## Columns: 5
## $ ID          <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ...
## $ name        <chr> "BRILLIANT TOTS PTE. LTD.", "BUBBLESLAND PLAYHOUSE PTE ...
## $ address     <chr> "610, JURONG WEST STREET 65, #01 - 534, S 640610", "238...
## $ postal_code <chr> "640610", "540238", "737856", "730369", "542327", "8212...
## $ geometry    <POINT [m]> POINT (13258.34 35611.04), POINT (35272.09 41373....

Import distance matrix

dm_preschools <- read_csv('data/aspatial/distance matrix/hdb_preschools.csv')

## Parsed with column specification:
## cols(
##   origin_id = col_double(),
##   destination_id = col_double(),
##   entry_cost = col_double(),
##   network_cost = col_double(),
##   exit_cost = col_double(),
##   total_cost = col_double()
## )

glimpse(dm_preschools)

## Rows: 18,876,550
## Columns: 6
## $ origin_id      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ destination_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ entry_cost     <dbl> NA, 89.29121, 89.29121, 89.29121, 89.29121, 89.29121...
## $ network_cost   <dbl> NA, 19008.163, 19731.005, 17167.844, 19267.844, 1926...
## $ exit_cost      <dbl> NA, 57.89667, 81.98225, 72.21746, 116.03101, 116.031...
## $ total_cost     <dbl> NA, 19155.351, 19902.279, 17329.353, 19473.166, 1947...

Tidy matrix – make wider

dm_preschools <- dm_preschools %>%
  # Obtain final matrix
  select(origin_id, destination_id, total_cost) %>%
  pivot_wider(names_from = destination_id, values_from = total_cost) %>%
  select(-origin_id)
  
# Convert to km
dm_preschools <- as.matrix(dm_preschools/1000)

Calculate average pre-school capacity

• Pre-schools include childcare centres and kindergartens.
• Average capacity for a pre-school is estimated using the total capacity of pre-schools in Singapore, divided by the total number of pre-schools in Singapore.
• Form the data, the total number of pre-schools in Singapore is 1925.
• Total capacity of pre-schools (childcare and kindergartens) is 215,579

preschool_capacity <- round(215579 / nrow(preschools))
preschool_capacity

## [1] 112

Calculate Hansen accessibility

preschool_capacity_list <- rep(preschool_capacity, nrow(preschools))

acc_Hansen <- data.frame(ac(hdb$children_each_hdb,
                  preschool_capacity_list,
                  dm_preschools, 
                  d0 = 50,
                  power = 2, 
                  family = "Hansen"))

colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hansen_preschools <- bind_cols(hdb, acc_Hansen)

hansen_preschools

hansen_preschools[is.na(hansen_preschools)] <- 0

glimpse(hansen_preschools)

## Rows: 9,806
## Columns: 11
## $ ID                <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15...
## $ house_number      <chr> "21", "21", "94B", "94B", "94C", "94C", "25B", "2...
## $ postal_code       <chr> "270021", "270021", "461094", "461094", "462094",...
## $ street            <chr> "Ghim Moh Road", "Ghim Moh Road", "Bedok North Av...
## $ num_levels        <chr> "10", "10", "13", "13", "13", "13", "27", "27", "...
## $ SUBZONE_N         <chr> "0", "0", "0", "0", "0", "0", "0", "0", "0", "0",...
## $ children_hdb      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ num_hdb           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ children_each_hdb <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ accHansen         <dbl> 0.00000, 316.59300, 36.35925, 36.56049, 191.26065...
## $ geometry          <POINT [m]> POINT (23015.11 32522.95), POINT (23015.11 ...

Function to compute risk weights

compute_risk_weights <- function(df, name) {
  first_quantile_hansen <- quantile(df$accHansen, 0.25)[[1]]
  median_hansen <- median(df$accHansen)
  third_quantile_hansen <- quantile(df$accHansen, 0.75)[[1]]
  max_hansen <- max(df$accHansen)

  df1 <- df %>%
    mutate('risk_{name}' := case_when(
      accHansen < first_quantile_hansen ~ 4,
      (accHansen >= first_quantile_hansen) & (accHansen < median_hansen) ~ 3,
      (accHansen >= median_hansen) & (accHansen < third_quantile_hansen) ~ 2,
      (accHansen >= third_quantile_hansen) ~ 1))
  
  return(df1)
}

Compute risk weights

hansen_cib_gardens <- compute_risk_weights(hansen_cib_gardens, 'cib_gardens')
hansen_dus_school_sports_facilities <- compute_risk_weights(hansen_dus_school_sports_facilities, 'dus_sports')
hansen_preschools <- compute_risk_weights(hansen_preschools, 'preschools')
hansen_pri_schools <- compute_risk_weights(hansen_pri_schools, 'pri_schools')
hansen_student_care <- compute_risk_weights(hansen_student_care, 'student_care')
hansen_watersports_facilities <- compute_risk_weights(hansen_watersports_facilities, 'watersports_facilities')

hdb_risk <- hdb %>%
  cbind(hansen_cib_gardens$risk_cib_gardens,
        hansen_dus_school_sports_facilities$risk_dus_sports,
        hansen_preschools$risk_preschools,
        hansen_pri_schools$risk_pri_schools,
        hansen_student_care$risk_student_care,
        hansen_watersports_facilities$risk_watersports_facilities) %>%
  rename(risk_cib_gardens = hansen_cib_gardens.risk_cib_gardens,
         risk_dus_sports = hansen_dus_school_sports_facilities.risk_dus_sports,
         risk_preschools = hansen_preschools.risk_preschools,
         risk_pri_schools = hansen_pri_schools.risk_pri_schools,
         risk_student_care = hansen_student_care.risk_student_care,
         risk_watersports_facilities = hansen_watersports_facilities.risk_watersports_facilities) %>%
  mutate(physical_health_risk = (risk_dus_sports + risk_watersports_facilities) / 2) %>%
  mutate(social_competence_risk = (risk_dus_sports + risk_preschools + risk_pri_schools + risk_student_care + risk_watersports_facilities) / 5) %>%
  mutate(emotional_maturity_risk = risk_cib_gardens) %>%
  mutate(composite_risk = (physical_health_risk + social_competence_risk + emotional_maturity_risk) / 3)

glimpse(hdb_risk)

## Rows: 9,806
## Columns: 20
## $ ID                          <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, ...
## $ house_number                <chr> "21", "21", "94B", "94B", "94C", "94C",...
## $ postal_code                 <chr> "270021", "270021", "461094", "461094",...
## $ street                      <chr> "Ghim Moh Road", "Ghim Moh Road", "Bedo...
## $ num_levels                  <chr> "10", "10", "13", "13", "13", "13", "27...
## $ SUBZONE_N                   <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ children_hdb                <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ num_hdb                     <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ children_each_hdb           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ risk_cib_gardens            <dbl> 3, 3, 3, 3, 2, 2, 1, 1, 2, 2, 2, 2, 2, ...
## $ risk_dus_sports             <dbl> 3, 3, 3, 3, 2, 2, 3, 3, 3, 3, 1, 1, 1, ...
## $ risk_preschools             <dbl> 4, 3, 4, 4, 4, 4, 3, 3, 2, 3, 3, 2, 2, ...
## $ risk_pri_schools            <dbl> 2, 2, 3, 3, 3, 3, 4, 4, 1, 2, 2, 1, 1, ...
## $ risk_student_care           <dbl> 3, 3, 2, 2, 3, 3, 3, 3, 2, 2, 2, 1, 1, ...
## $ risk_watersports_facilities <dbl> 1, 1, 2, 2, 2, 2, 4, 4, 1, 1, 1, 1, 1, ...
## $ geometry                    <POINT [m]> POINT (23015.11 32522.95), POINT ...
## $ physical_health_risk        <dbl> 2.0, 2.0, 2.5, 2.5, 2.0, 2.0, 3.5, 3.5,...
## $ social_competence_risk      <dbl> 2.6, 2.4, 2.8, 2.8, 2.8, 2.8, 3.4, 3.4,...
## $ emotional_maturity_risk     <dbl> 3, 3, 3, 3, 2, 2, 1, 1, 2, 2, 2, 2, 2, ...
## $ composite_risk              <dbl> 2.533333, 2.466667, 2.766667, 2.766667,...

8 Visualise risk score

tm_shape(subzone_children) +
  tm_polygons(col = 'gray90') +
  tm_shape(hdb_risk) +
  tm_bubbles(col = 'composite_risk', 
             size = 0.1,
             alpha = 0.5,
             border.lwd = NA)

9 Save data files to shiny app

# st_write(hdb_risk, '../shinyapp/data/geospatial/hdb_risk.shp', driver = 'ESRI Shapefile')

# st_write(subzone_children, '../shinyapp/data/geospatial/subzone_children.shp', driver = 'ESRI Shapefile')

tmap_mode('view')

## tmap mode set to interactive viewing

tm_shape(hdb_risk) +
  tm_bubbles()

tmap_mode('plot')

## tmap mode set to plotting

hdb_accessibility <- hdb %>%
  cbind(hansen_cib_gardens$accHansen,
        hansen_dus_school_sports_facilities$accHansen,
        hansen_preschools$accHansen,
        hansen_pri_schools$accHansen,
        hansen_student_care$accHansen,
        hansen_watersports_facilities$accHansen) %>%
  rename(hansen_cib_gardens = hansen_cib_gardens.accHansen,
         hansen_dus_sports = hansen_dus_school_sports_facilities.accHansen,
         hansen_preschools = hansen_preschools.accHansen,
         hansen_pri_schools = hansen_pri_schools.accHansen,
         hansen_student_care = hansen_student_care.accHansen,
         hansen_watersports_facilities = hansen_watersports_facilities.accHansen)
hdb_accessibility

10 Export accessibility data file

# st_write(hdb_accessibility, '../shinyapp/data/geospatial/hdb_accessibility.shp', driver = 'ESRI Shapefile')

tmap_mode('view')

## tmap mode set to interactive viewing

tm_shape(hdb_accessibility) +
  tm_bubbles(col = 'hansen_student_care')